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We perform numerical simulations of purely repulsive soft colloidal particles interacting via a 
generalized elastic potential and constrained to a two-dimensional plane and to the surface of a 
spherical shell. For the planar case, we compute the phase diagram in terms of the system's rescaled 
density and temperature. We find that a large number of ordered phases becomes accessible at low 
temperatures as the density of the system increases, and we study systematically how structural 
variety depends on the functional shape of the pair potential. For the spherical case, we revisit the 
generalized Thomson problem for small numbers of particles A'^ < 12 and identify, enumerate and 
compare the minimal energy polyhedra established by the location of the particles to those of the 
corresponding electrostatic system. 



Introduction 

Understanding how nanocomponents spontaneously 
organize into complex macroscopic structures is one of 
the great challenges in the field of soft matter today. In- 
deed, the ability to predict and control the phase behav- 
ior of a solution, given a set of components, may open the 
way to the development of materials with novel optical, 
mechanical, and electronic properties. Despite much ef- 
fort in this direction, the question of how to link the phys- 
ical character of the components, i.e. their bare shape 
and their interaction potential, to their phase behavior 
still remains unanswered. 

Of particular interest for optical applications are the 
crystalline phases. Until recently, it was believed that 
shape anisotropy and/or directional interactions where 
key elements for the formation of crystals with complex 
symmetries commonly observed for instance in atomic 
solids. The study of the phase behavior of components in- 
teracting with exotic (yet isotropic) potentials has proven 
this belief to be incorrect. For instance, Torquato et al. 
have shown, using inverse optimization techniques, that 
it is possible to achieve non-closed-packed structures such 
as simple cubic, hexagonal, wurtzite and even diamond 
phases with isotropic pair potentials (see ref for a re- 
view on the subject). Furthermore, spherical particles 
interacting via a hard core and a repulsive shoulder po- 
tential are known to organize into complex mesophases 
not unlike those found with diblock copolymers 

One class of pair potentials that has recently attracted 
much attention is that describing the interactions be- 
tween soft/deformable mesoparticles. Unlike typical col- 
loidal particles for which excluded volume interactions 
are strictly enforced via a hard-core or a Lennard-Jones 
potential, complex mesoparticles such as charged or neu- 
tral star polymers, dendrimers or microgels present a 
more peculiar pair potential describing their volume in- 
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teractions. Surprisingly, the simple relaxation of the 
constraint of mutual impenetrability between isotropic 
components gives access to several non-close-packed crys- 
talline structures at high densities. Given the complex- 
ity of these mesoparticles, their interactions are usually 
extracted via an explicit coarse-grained procedure to ob- 
tain ad-hoc effective pair potentials. What emerges is 
a great variety of nontrivial interactions, some of which 
allow for even complete overlap among the components, 
resulting in a very rich phenomenological behavior. Re- 
markably, it is feasible to engineer interactions between 
star polymers or dendrimers by controlling their overall 
chemical/topological properties 

The phase behavior of several systems adopting these 
exotic, but physically inspired, interactions has been the 
subject of several publications [4Hl4]. Notably, it was 
found that some classes of soft interactions lead to reen- 
trant melting transitions, others to polymorphic clus- 
ter phases |15j, and in general to multiple transitions 
involving close-packed and non-close-packed crystalline 
phases [TMT5| as a function of the system density. Re- 
markably, the phase behavior of these systems is very 
much dependent on the the shape of the pair poten- 
tial. Likos et al. [19] established a criterion to predict 
whether for a bounded and repulsive potential reentrant 
melting or cluster phases will occur based on the sign of 
the Fourier transform of the interaction. Nevertheless, 
there is currently no method to predict a priori what 
specific crystal structures may become accessible for a 
given potential. For a recent review on the subject we 
refer the reader to reference [20] . 

In this paper we focus on two-dimensional systems. 
Although much of the research in this field has focused 
on three-dimensional systems, earlier numerical work on 
particles interacting via a hard core plus a shoulder po- 
tential in two dimensions has revealed a phenomenologi- 
cal behavior that can be as rich as that observed at higher 
dimensionality pT1124) . Here we focus on bounded soft 
potentials. Specifically, expanding on our recent results 
on Hertzian spheres [18] and dumbbells [25] in three di- 
mensions, we analyze the phase behavior of soft elastic 
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particles in two dimensions. We generalize our results 
to spherical surfaces, and discuss the generalized Thom- 
son problem for our soft potentials by identifying novel 
polyhedral structures representing minimal energy con- 
figurations formed by soft particles constrained on the 
surface of a sphere at different packing densities. 

Methods 

To study the phase behavior of soft nanoparticles in 
two dimensions we performed numerical simulations. We 
considered systems of at least N = 1000 spherical parti- 
cles and used the standard Monte Carlo method in the 
NPT and NVT ensemble for the planar and spherical case 
respectively. All simulations where run for a minimum 
of 10^ iterations. Any two particles in our system inter- 
act via a Hertz potential. The Hertz potential describes 
the elastic energy penalty associated with an axial com- 
pression of two deformable spheres. Its functional form 
can be generalized as follows 

I lor r > a 

where cr is the particle diameter, r is the interparticle 
distance, e is the unit of energy, and a is a parameter 
that we control to modulate the shape of the interaction. 
The elastic case is recovered for a ~ 5/2. 

The Hertz model, developed to account for small elas- 
tic deformations, becomes inaccurate when associated to 
the large overlaps among particles that is achieved at 
large densities; nevertheless, it provides us with a simple 
representation of a finite-ranged, bounded soft potential 
with a positive definite Fourier transform (a condition 
that guarantees reentrant melting in the phase diagram.) 
We included in our study two more values of a to account 
for a slightly harder {a = 3/2) and a slightly weaker 
{a = 7/2) potential than the elastic one (a = 5/2.) Fig- 
ure [T] shows these potentials. 

Crystallization in our simulations was determined us- 
ing a combination of visual inspection and numerical or- 
der parameters similar to those described in . Specifi- 
cally, given a randomly chosen particle j , and a reference 
direction set to be along the bond between this parti- 
cle and one of its nearest neighbors, we can define an 
order parameter sensitive to a crystal of D-fold symme- 
try (specifically six-fold hexagonal and four-fold square 
crystals) as follows: For all particles k within a radius of 
some cutoff Vc of j (including j itself), we calculate the 
angle (pik between the bond joining each of the D nearest 
neighbors / of A: to fc and the reference direction. If the 
number of particles within of j is N^. , then the order 
parameter Vl/j^ is given by: This process is then repeated 
N times and an average ^'^i is calculated, yielding a value 
between (for completely disordered systems) and 1 (for 
perfectly crystalline systems). This order parameter is 
analogous to that defined in [Mj except that instead of 
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FIG. 1: Plot of the rescaled pair potential 
V{r)/e = (1 — ^)" for the three values of a used in this 
study: a= {f,|,|} 

being an average over all particles k in the system, it 
is only over particles k within Vc (= 1.65f7 in our case); 
thus it is "semi-local" and allows crystalline order to be 
detected even in the presence of metastable grain bound- 
aries. The order parameter in |26j represents the rc — > oo 
limit of the above method. 

To draw the phase diagrams we scan the phase space 
defined by the rescaled temperature = k^T/e and the 
rescaled number density p,- = Na'^/A (where A is the 
area of the system) to evaluate the several phases of the 
system in terms (when possible) of its order parameter 
^ D . The lines separating the different regions are guides 
to the eye. Given the peculiarities of phase transitions in 
two dimensions, and the richness in the phases behavior 
reported, we have not attempted to establish equilibrium 
boundaries between the different phases, yet we have re- 
produced our results using Molecular Dynamics simula- 
tions with a Langevin Thermostat performed in key re- 
gions of the phase space, and compared the stability of 
the crystalline structures by measuring their relative en- 
ergies at Tr 0. 



Results 

We begin by presenting the phase diagram for the case 
of Hertzian (elastic) spheres. The result is shown in 
Fig. [2] The labels in the fi gure refer to the structures 
shown in Fig. |4] As expected from Likos's criterion [T9] 
our system shows reentrant melting behavior and a high- 
est temperature above which the system remains fluid at 
all densities. 

The phase behavior here reported is qualitatively sim- 
ilar to that observed in three dimensions, however the 
two-dimensional system does not produce a phase space 
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FIG. 2: Temperature vs. number density phase 
diagrams for two-dimensional particles interacting via a 
Hertz potential (a — 5/2). Labels on the plots refer to 
the labels of the phases in Fig. |4j 



that is as structurally rich as its three dimensional coun- 
terpart. Across the spectrum of densities (up to p,- = 
10) and temperatures (down to Tj- = 10~^) that we 
explored, we find multiple crystal-to-crystal transitions 
from hexagonal (phase a) to square symmetry (phase b) . 
Interestingly, the occurrence of the two phases has an 
alternating periodic pattern, and no isostructural transi- 
tion was detected in our system [27] . 

This pattern is only broken at a small range of densi- 
ties centered around pr — 3.25. This region (phase h) is 
characterized by a complex structure dominated by par- 
ticles arranged mostly into pentagonal units. Our color 
coding of the phase shown in Fig. |4]jh) helps rationaliz- 
ing the overall symmetry of the structure which in this 
representation can be thought of as the combination of 
four square lattices. Three lattices are generated by pen- 
tagonal units, displaced with respect to each other but 
oriented along the same axis. The fourth lattice is built 
out of the non-connected particles and is oriented along 
an axis that forms an angle of 45 degrees with the other 
lattices. 

To understand how the phase behavior is affected by 
the specific choice of the functional form of the poten- 
tial, we repeated the calculation of the phase diagram for 
a = 7/2 and a = 3/2. The results are presented is Fig.js] 
Although the overall trend is very similar, the number 
and the sequence of phases that we find are significantly 
different. Namely, a = 3/2 leads to a significantly richer 
structural behavior, while the large value of a results in a 
single hexagonal crystalline structure. As a reference, no- 
tice that our potential tends to the linear-ramp potential 
for a = 1 (for which a large number of phases including 
quasi-crystals have been observed ^3\) , we recover the 
square-shoulder potential for a — >• (for which a cascade 
of cluster phases are expected) , and tends to a Kronecker 
delta potential for a — > cx) (for which no crystallization 
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FIG. 3: Temperature vs. number density phase 
diagrams for two-dimensional particles interacting via 
the generalized Hertz potential with a — \ (top) and 
a = I (bottom). Labels on the plots refer to the labels 
of the phases in Fig. |4[ 



is expected at any finite density). 

At least as a general trend, it is therefore not surpris- 
ing that structural variety increases for smaller values of 
a. In this regard, it is interesting to estimate the criti- 
cal power above which the hexagonal lattice becomes the 
only stable lattice as observed in Fig. [sf^bottom). Given 
that the only other competing structure at large values 
of a is the square lattice, we computed the energy den- 
sity difference between hexagonal and square lattices as a 
function of density and at zero temperature for different 
values of a. The result is presented in Fig. [5] for densi- 
ties up to Pr = 25 (our full numerical data extend up to 
Pt- = 100, but are not shown for the sake of clarity). We 
find that the hexagonal lattice becomes more stable than 
the square lattice at all densities for any powers above the 
onset value a* ~ 3.182. Several numerical simulations at 
finite densities for powers ranging from a = 10 to a = 50 
where carried out, and indeed only hexagonal structures 



(a) Hexagonal lattice 



(b) Square lattice 



(c) Lines 



(d) Stretched honeycomb 
lattice 




(e) Kagome lattice (f) Open square lattice (g) Snub square tiling (h) Pentagons 

FIG. 4: Two-dimensional crystal phases formed in these systems. 
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FIG. 5: Energy density difference between hexagonal 
and square lattice at Tj. = as a function of system 
density at different values of a. We estimated the 
critical a to be around a ~ 3.1822. 



where observed upon system ordering. The presence of 
reentrant melting was also established for a = 5 and 
a = 10. The oscillating behavior of the energy density 
difference for a < a* also explains qualitatively the pe- 
riodic structural pattern observed in the phase diagram 
for a = 5/2. 

We now turn to the discussion of the second part 
of the paper. Given the structural variety observed in 
the two dimensional compression of soft nanoparticles, 



it is interesting to extend our results to other geome- 
tries. The spherical case is of particular interest be- 
cause a large body of work has been dedicated to find- 
ing and enumerating minimal energy conformations of 
particles constrained over a spherical shell, and repelling 
each other via an electrostatic repulsion [28l - [39| . This 
is what Thomson asked in 1904 [40] when attempting 
to construct his plum pudding model of the atom. Al- 
though there is consensus on the minimal energy struc- 
tures for TV < 100, the problem becomes more involved 
for large values of N for which an exponentially large 
number of low energy configurations becomes available. 
In this regime, lattices with overall icosahedral symmetry 
(icosadeltahedra) have been initially postulated as possi- 
ble global minima for specific magic number of parti- 
cles [Mj, but eventually it was shown that it is possible 
to lower the energy of these configurations adding dislo- 
cation defects [37l |4T| . Grain boundary scars have also 
been observed in experiments |42| . 

Although we confirm that the same pattern of phases 
obtained in two dimensions also develops on large areas 
of spherical shells - apart from the topologically required 
defects - as a function of the radius of the sphere (see 
Fig. [g] for a sample of these configurations) , we have not 
attempted to enumerate all possible low energy states 
for large N. It is however clear that considering soft po- 
tentials introduces two new parameters to the system: 
namely the packing density and the shape of the poten- 
tial. A careful analysis of the defects and symmetries 
arising in this regime is currently under investigation and 
will be the subject of another publication; here we focus 
on the case of small number of particles iV < 12 for 



5 




a = 3/2 and a = 5/2. 

Figures [7] and |8] show the results of our analysis for 
N G {5, 12}. The rescaled configurational energy E de- 
fined as the total internal energy divided by e, the re- 
duced radius tq = R/cr of the spherical surface and the 
names of the polyhedra (when available) are also given. 
The reference configurations for the electrostatic problem 
are highlighted with a dark frame. These data are ob- 
tained using both Monte Carlo and Molecular Dynamics 
simulations with a standard temperature annealing pro- 
cedure down to Tr — )• 0. 

Unlike the case of charged point particles, we find that 
indeed the number and the nature of the structures is 
very much dependent on the radius of the sphere. The 
overall trend follows what found in the analysis of the pla- 
nar two-dimensional system; namely, smaller values of a 
lead to richer structural diversity. In all cases we where 
able to recover the global minima of the electrostatic sys- 
tem for at least one spherical radius. Interestingly, = 4 
and iV = 6 are the only cases in which for both powers 
we obtain only one structure, a tetrahedron and a octa- 
hedron respectively (not shown), which are stable for any 
value of r-Q explored (down to tq = 0.10). The stability 
of the octahedron {N = 6) was also confirmed for values 
of a up to 10. Low energy configurations for values of N 
larger than 12 have also been considered, however, given 
the large variety and the complexity of the structures 
arising upon increasing iV, we have not attempted to 
enumerate them. In fact, more sophisticated algorithms 
than the one used in this work would be necessary to 
thoroughly explore the energy landscape in search of a 
global minimum. 

Conclusions 

In this paper we explore the phase behavior of col- 
loidal particles interacting via a generalized soft (Hertz) 
potential. Specifically, we analyze how structural diver- 
sity depends on the particular choice of the functional 
form of the potential for particles constrained on a two- 



dimensional plane and on the surface of a spherical shell. 
For the planar case we compute how the phase diagram 
of the system changes with the functional form of the 
pair potential, and establish a limit above which (for the 
class of potentials explored) hexagonal packing becomes 
the only allowed symmetry of the ordered phase. For 
small number of particles [N < 12), we identify on a 
spherical shell the polyhedral configurations establishing 
global energy minima for different pair potentials, and 
compare the results to the global minima corresponding 
to the classic Thomson problem. We show that unlike 
the electrostatic case, for the same number of particles 
the system presents more than one global minimum de- 
pending on the radius of the sphere. 

It would be interesting to study the geometry of the 
defects arising for larger number of particles on the spher- 
ical shell, analyze the overall symmetry of the ground 
states, and explore whether patterns of grain boundaries, 
expected for large values of N in the electrostatic sys- 
tem, would also manifest for soft potentials. This is a 
very challenging problem that we have already begun to 
investigate and requires more sophisticated minimization 
techniques than the ones employed in this study. For the 
time being, we have established that the same phases ob- 
served in the planar geometry do develop on the spherical 
shell at comparable number densities. 

Finally, given the sensitivity of the phase behavior on 
the specific choice of the pair potential, it becomes critical 
at this point to develop sophisticated and reliable coarse- 
graining procedures able to capture accurate interactions 
between dendrimers or star polymers, but also the inverse 
problem, i.e. the design of polymeric structures able to 
reproduce desired functional forms for the pair potential, 
is in great need of a better theoretical understanding. 
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Cluster geometry 





square pyramid 
ro = 0.65 
E = 0.116675 



p(6) 

ro = 0.58 
E = 0.457195 




ro = 0.75 
E = 0.1481239 




square antiprism 
ro = 0.70 
E = 0.851626 



triaugmented 
triangular prism 
ro = 0.81 
E = 0.268064 




pentagonal 
dipyramid 
ro = 0.65 
E = 0.802319 




sphenocorona 

ro = 0.90 
E = 0.035971 




gyroclongatcd 
pentagonal pyramid 
ro = 0.90 
E = 0.244728 




icosahedron 

ro = 0.90 
E = 0.373165 




cube 
ro = 0.58 
E = 2.423538 




gyroelongated 
square pyramid 
ro = 0.78 
E = 0.508830 




pentagonal prism 

ro = 0.60 
E = 4.102489 




3(11) 



ro = 0.80 
E = 1.308288 




elongated 
pentagonal 
dipyramid 
ro = 0.75 
E = 2.857171 




trigonal dipyramid 
ro = 0.50 
E = 1.098194 







ro = 0.58 
E = 1.635489 



ro = 0.55 
E = 2.066014 



ro = 0.50 
E = 2.956885 



ro = 0.40 
E = 5.463620 




square antiprism 
ro = 0.50 
E = 4.212084 




triaugmented 
triangular prism 
ro = 0.70 
E = 1.424194 




pentagonal 
antiprism 
ro = 0.57 
E = 4.921405 




elongated square 
pyramid 
ro = 0.63 
E = 2.489422 




ro = 0.58 
E = 3.479048 




triaugmented 
triangular prism 
ro = 0.50 
E = 5.696124 




gyroelongated 
square dipyramid 
ro = 0.40 
E = 12.689921 






elongated 
pentagonal pyramid 
ro = 0.68 
E = 3.322424 



ro = 0.63 
E = 4.479854 



gyrclongatcd 
pentagonal pyramid 
ro = 0.55 
E = 7.042131 





p(12) 

ro = 0.70 
E = 3.806734 



p(12) 

ro = 0.65 
E = 5.087802 




icosahedron 

ro = 0.50 
E = 11.529991 




ro = 0.40 
E = 15.783736 



FIG. 7: Global minimum configurations formed by small number of particles constrained on the surface of a sphere 
interacting with a generalized hertzian soft potential for a = 3/2, at different values of reduced sphere radius ro- 
The reduced energy of the configurations is also indicated together with the name of the polyhedric structures. 
Some structures for which no name was found are indicated with the nomenclature ■ For A'^ = 4 and iV = 6 
(not shown) we obtain respectively a tetrahedral and an octahedral structure for every rg. Reference configurations 
for long-range potentials have been highlighted with a dark border. 
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Cluster geometry 




square pyramid 
ro = 0.70 
E = 0.000047 




ro = 0.75 
E = 0.008832 




triaugmented 
triangular prism 
ro = 0.75 
E = 0.110712 




ro = 0.58 
E = 0.081254 




pentagonal 
dipyramid 
ro = 0.50 
E = 1.012408 




gyroelongated 
square pyramid 
ro = 0.72 
E = 0.196371 




trigonal dipyramid 
ro = 0.50 
E = 0.298277 




ro = 0.25 
E = 6.627468 




triaugmented 
triangular prism 
ro = 0.50 
E = 2.203406 







sphcnocorona 

ro = 0.80 
E = 0.097219 



gyroelongated 

square dipyramid 
ro = 0.65 
E = 0.850254 



sphcnocorona 

ro = 0.60 
E = 1.361079 



gyroelongated 
square dipyramid 
ro = 0.50 
E = 2.986067 





gyroelongated 
pentagonal pyramid 
ro = 0.80 
E = 0.221753 



ro = 0.70 
E = 0.793120 



FIG. 8: Global minima configurations formed by small number of particles constrained on the surface of a sphere 
interacting with a generalized hertzian soft potential for a = 5/2, at different values of reduced sphere radius ro- 
The reduced energy of the configurations is also indicated together with the name of the polyhedric structures. 
Some structures for which no name was found are indicated with the nomenclature . For N = 4, N = 6 and 
A'' = 12 (not shown) we obtain respectively a tetrahedral, an octahedral and an icosahedral structure for every ro. 
Reference configurations for long-range potentials have been highlighted with a dark border. 
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